Let’s start by reading your data:
# Import libraries:
# Load libraries
# Note: commented-out libraries are experimental libraries for other functions not used in this final analysis file
#load some useful libraries
library(limma)
library(edgeR)
library(sns)
library(data.table)
library(tidyr)
library(tidyverse)
library(sjmisc)
require(foreign)
library(foreign)
require(ggplot2)
require(MASS)
require(Hmisc)
require(reshape2)
library(stringr)
library(rlist)
library(rlang)
library(forcats)
library(dplyr)
library(tidyverse)
library(ggh4x)
require(tidyselect)
library(DataCombine)
library(ggrepel)
library(ggbeeswarm)
library(ggsignif)
library(ggpubr)
library(btools)
library(ggpmisc)
library(quantreg)
library(broom)
library(gginnards)
library(broom.mixed)
library(scales)
#Corncob and Phyloseq:
library(corncob)
library(compositions)
Welcome to compositions, a package for compositional data analysis.
Find an intro with "? compositions"
Attaching package: ‘compositions’
The following object is masked from ‘package:SparseM’:
norm
The following objects are masked from ‘package:stats’:
anova, cor, cov, dist, var
The following objects are masked from ‘package:base’:
%*%, norm, scale, scale.default
library(phyloseq)
# Import data:
otus <- read.csv("gut_full.csv", check.names=F)
taxa <- read.csv("taxa.csv",sep='\t')
otus$sex <- factor(otus$sex) #factorize sex
otus$vendor_dashboard <- factor(otus$vendor_dashboard)
otus$bowel <- factor(as.numeric(factor(otus$bowel, levels = c(1,2,3,4), labels = c(1,2,3,4))))
otus <- otus[ which(otus$vendor_dashboard == "Second Genome" | otus$vendor_dashboard == "research-microbiome"),] # keep only data where vendor is explicit
otus$vendor_dashboard = str_replace_all(otus$vendor_dashboard,"research-microbiome","DNA Genotek")
otus$vendor_dashboard <- factor(otus$vendor_dashboard) # factorize vendor
#otus[ which(otus$vendor_dashboard == "Second Genome"),]
otus <- within(otus, bowel <- relevel(bowel, ref = 3))
# Pre-processing and checking validity of data
# Algorithms provided by Christian Diener, PhD:
##############################################################################
dat <- otus
otus <- otus[!duplicated(otus$public_client_id), ]
genus_cols <- grepl("taxa_", names(otus))
rownames(otus) <- otus$public_client_id
sdata <- otus[, !genus_cols]
bowel <- paste0(sdata['bowel'])
table(sdata['bowel'])
3 1 2 4
1796 83 769 41
otus <- as.matrix(otus[, genus_cols])
colnames(otus) <- gsub("taxa_", "", colnames(otus))
tax_matrix <- as.matrix(taxa[, 2:ncol(taxa)])
rownames(tax_matrix) <- taxa[, 1]
as.data.frame(otus)
Now we convert it to fit in a phyloseq object. We need to fulfill the following rules:
To double check if everything works let’s do some validity checks:
stopifnot(all(rownames(otus) %in% rownames(sdata)))
stopifnot(nrow(otus) == nrow(sdata))
print("sample names match")
[1] "sample names match"
stopifnot(all(colnames(otus) %in% rownames(tax_matrix)))
stopifnot(ncol(otus) == nrow(tax_matrix))
stopifnot(!anyDuplicated(tax_matrix))
print("taxa look okay")
[1] "taxa look okay"
stopifnot(!anyDuplicated(sdata))
print("sample data looks okay")
[1] "sample data looks okay"
If that passes we can go ahead and build our phyloseq object.
ps <- phyloseq(
otu_table(otus, taxa_are_rows = FALSE),
tax_table(tax_matrix),
sample_data(sdata)
)
ps
phyloseq-class experiment-level object
otu_table() OTU Table: [ 68 taxa and 2689 samples ]
sample_data() Sample Data: [ 2689 samples by 7 sample variables ]
tax_table() Taxonomy Table: [ 68 taxa by 6 taxonomic ranks ]
names(sample_data(ps))
[1] "public_client_id" "sex" "age" "BMI_CALC"
[5] "vendor_dashboard" "bowel" "eGFR"
bowel <- sample_data(ps)$bowel
##############################################################################
#Differential Test
#Uncomment this code chunk to run the CORNCOB modeling.
#Otherwise, load the dv_analysis object already saved from a previously completed modeling run.
#dv_analysis <- differentialTest(formula = ~ bowel + sex + age + BMI_CALC + vendor_dashboard + eGFR,
# phi.formula = ~ 1,
# formula_null = ~ sex + age + BMI_CALC + vendor_dashboard + eGFR,
# phi.formula_null = ~ 1,
# data = ps,
# test = "LRT", boot = FALSE,
# full_output = TRUE,
# fdr_cutoff = 0.05)
#Load a previous CORNCOB modeling
#Click on the "corncob.rds" file in the folder with R Studio open and add the dataframe to the environment. Name it "corncob":
dv_analysis <- corncob
#See the significant taxa from the computation:
dv_analysis$significant_taxa
[1] "Bacteria.Actinobacteria.Actinobacteria.Bifidobacteriales.Bifidobacteriaceae.Bifidobacterium"
[2] "Bacteria.Actinobacteria.Coriobacteriia.Coriobacteriales.Coriobacteriaceae.Collinsella"
[3] "Bacteria.Actinobacteria.Coriobacteriia.Coriobacteriales.Eggerthellaceae.Adlercreutzia"
[4] "Bacteria.Bacteroidetes.Bacteroidia.Bacteroidales.Bacteroidaceae.Bacteroides"
[5] "Bacteria.Bacteroidetes.Bacteroidia.Bacteroidales.Marinifilaceae.Odoribacter"
[6] "Bacteria.Bacteroidetes.Bacteroidia.Bacteroidales.Rikenellaceae.Alistipes"
[7] "Bacteria.Bacteroidetes.Bacteroidia.Bacteroidales.Tannerellaceae.Parabacteroides"
[8] "Bacteria.Firmicutes.Bacilli.Lactobacillales.Streptococcaceae.Streptococcus"
[9] "Bacteria.Firmicutes.Clostridia.Clostridiales.Christensenellaceae.Christensenellaceae_R-7_group"
[10] "Bacteria.Firmicutes.Clostridia.Clostridiales.Christensenellaceae.nan"
[11] "Bacteria.Firmicutes.Clostridia.Clostridiales.Clostridiales_vadinBB60_group.nan"
[12] "Bacteria.Firmicutes.Clostridia.Clostridiales.Family_XIII.Family_XIII_AD3011_group"
[13] "Bacteria.Firmicutes.Clostridia.Clostridiales.Family_XIII.Family_XIII_UCG-001"
[14] "Bacteria.Firmicutes.Clostridia.Clostridiales.Family_XIII.nan"
[15] "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.Agathobacter"
[16] "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.Blautia"
[17] "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.Lachnoclostridium"
[18] "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.Lachnospira"
[19] "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.Lachnospiraceae_NK4A136_group"
[20] "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.Lachnospiraceae_UCG-004"
[21] "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.Marvinbryantia"
[22] "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.nan"
[23] "Bacteria.Firmicutes.Clostridia.Clostridiales.Peptostreptococcaceae.Romboutsia"
[24] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Anaerotruncus"
[25] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Butyricicoccus"
[26] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Candidatus_Soleaferrea"
[27] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.DTU089"
[28] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Faecalibacterium"
[29] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.GCA-900066225"
[30] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Intestinimonas"
[31] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Negativibacillus"
[32] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Oscillibacter"
[33] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminiclostridium_5"
[34] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminiclostridium_9"
[35] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminococcaceae_NK4A214_group"
[36] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminococcaceae_UCG-002"
[37] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminococcaceae_UCG-004"
[38] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminococcaceae_UCG-005"
[39] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminococcaceae_UCG-013"
[40] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminococcus_1"
[41] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminococcus_2"
[42] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Subdoligranulum"
[43] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.UBA1819"
[44] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.nan"
[45] "Bacteria.Firmicutes.Erysipelotrichia.Erysipelotrichales.Erysipelotrichaceae.Holdemania"
[46] "Bacteria.Firmicutes.Erysipelotrichia.Erysipelotrichales.Erysipelotrichaceae.nan"
[47] "Bacteria.Verrucomicrobia.Verrucomicrobiae.Verrucomicrobiales.Akkermansiaceae.Akkermansia"
#prepare the taxa df for manipulation:
taxa<-dv_analysis$all_models
dv_analysis$significant_taxa
[1] "Bacteria.Actinobacteria.Actinobacteria.Bifidobacteriales.Bifidobacteriaceae.Bifidobacterium"
[2] "Bacteria.Actinobacteria.Coriobacteriia.Coriobacteriales.Coriobacteriaceae.Collinsella"
[3] "Bacteria.Actinobacteria.Coriobacteriia.Coriobacteriales.Eggerthellaceae.Adlercreutzia"
[4] "Bacteria.Bacteroidetes.Bacteroidia.Bacteroidales.Bacteroidaceae.Bacteroides"
[5] "Bacteria.Bacteroidetes.Bacteroidia.Bacteroidales.Marinifilaceae.Odoribacter"
[6] "Bacteria.Bacteroidetes.Bacteroidia.Bacteroidales.Rikenellaceae.Alistipes"
[7] "Bacteria.Bacteroidetes.Bacteroidia.Bacteroidales.Tannerellaceae.Parabacteroides"
[8] "Bacteria.Firmicutes.Bacilli.Lactobacillales.Streptococcaceae.Streptococcus"
[9] "Bacteria.Firmicutes.Clostridia.Clostridiales.Christensenellaceae.Christensenellaceae_R-7_group"
[10] "Bacteria.Firmicutes.Clostridia.Clostridiales.Christensenellaceae.nan"
[11] "Bacteria.Firmicutes.Clostridia.Clostridiales.Clostridiales_vadinBB60_group.nan"
[12] "Bacteria.Firmicutes.Clostridia.Clostridiales.Family_XIII.Family_XIII_AD3011_group"
[13] "Bacteria.Firmicutes.Clostridia.Clostridiales.Family_XIII.Family_XIII_UCG-001"
[14] "Bacteria.Firmicutes.Clostridia.Clostridiales.Family_XIII.nan"
[15] "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.Agathobacter"
[16] "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.Blautia"
[17] "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.Lachnoclostridium"
[18] "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.Lachnospira"
[19] "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.Lachnospiraceae_NK4A136_group"
[20] "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.Lachnospiraceae_UCG-004"
[21] "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.Marvinbryantia"
[22] "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.nan"
[23] "Bacteria.Firmicutes.Clostridia.Clostridiales.Peptostreptococcaceae.Romboutsia"
[24] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Anaerotruncus"
[25] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Butyricicoccus"
[26] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Candidatus_Soleaferrea"
[27] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.DTU089"
[28] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Faecalibacterium"
[29] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.GCA-900066225"
[30] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Intestinimonas"
[31] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Negativibacillus"
[32] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Oscillibacter"
[33] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminiclostridium_5"
[34] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminiclostridium_9"
[35] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminococcaceae_NK4A214_group"
[36] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminococcaceae_UCG-002"
[37] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminococcaceae_UCG-004"
[38] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminococcaceae_UCG-005"
[39] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminococcaceae_UCG-013"
[40] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminococcus_1"
[41] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminococcus_2"
[42] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Subdoligranulum"
[43] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.UBA1819"
[44] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.nan"
[45] "Bacteria.Firmicutes.Erysipelotrichia.Erysipelotrichales.Erysipelotrichaceae.Holdemania"
[46] "Bacteria.Firmicutes.Erysipelotrichia.Erysipelotrichales.Erysipelotrichaceae.nan"
[47] "Bacteria.Verrucomicrobia.Verrucomicrobiae.Verrucomicrobiales.Akkermansiaceae.Akkermansia"
as.data.frame(dv_analysis$p_fdr)
taxa_p <- DropNA(as.data.frame(dv_analysis$p_fdr), Var = "dv_analysis$p_fdr")
0 rows dropped from the data frame because of missing values.
taxa_p[1]
[1] 0.02120539
taxa[[1]]$coefficients
Estimate Std. Error t value Pr(>|t|)
mu.(Intercept) -3.368074680 0.182186294 -18.4869816 6.486060e-72
mu.bowel1 0.108766498 0.121122817 0.8979852 3.692741e-01
mu.bowel2 0.136368966 0.046755183 2.9166599 3.567424e-03
mu.bowel4 -0.212360465 0.178644229 -1.1887340 2.346497e-01
mu.sexW -0.212789194 0.048369795 -4.3992163 1.128862e-05
mu.age -0.016148200 0.001821779 -8.8639749 1.378229e-18
mu.BMI_CALC -0.003154924 0.003370311 -0.9360928 3.493098e-01
mu.vendor_dashboardSecond Genome 0.127694322 0.060993339 2.0935782 3.639137e-02
mu.eGFR 0.001582766 0.001243639 1.2726892 2.032388e-01
phi.(Intercept) -2.808912916 0.043294681 -64.8789385 0.000000e+00
#Load the arivale_phylo.rds file in the project folder and load the object into the environment.
#Rarefy the genotek dataset to an even depth using phyloseq.
sample_data(arivale_phylo)
Sample Data: [3653 samples by 26 sample variables]:
#Tom Wilmanski, PhD's code:
########################################################################################
rarefied_genotek=rarefy_even_depth(arivale_phylo, sample.size = min(sample_sums(arivale_phylo)),
rngseed = 111, replace = FALSE, trimOTUs = TRUE, verbose = TRUE)
`set.seed(111)` was used to initialize repeatable random subsampling.
Please record this for your records so others can reproduce.
Try `set.seed(111); .Random.seed` for the full vector
...
25984OTUs were removed because they are no longer
present in any sample after random subsampling
...
richness <- estimate_richness(rarefied_genotek, measures=c("Shannon","Observed"))
saveRDS(richness, "richness.rds")
########################################################################################
richness$public_client_id <-sample_data(rarefied_genotek)$public_client_id
library(btools)
richness$Pielou <- richness$Shannon/richness$Observed
richness
hist(richness$Pielou,
main = "Pielou's Evenness (Shannon/Observed ASV)",
xlab = "Evenness",
xlim = c(0,0.05),
breaks = 200)
hist(richness$Shannon,
main = "Shannon Diversity",
xlab = "Diversity Index",
xlim = c(0,6),
breaks = 20)
hist(richness$Observed,
main = "Observed ASVs",
xlab = "ASVs",
xlim = c(0,1000),
breaks = 30)
combinations <- list(c("Constipation","Low Normal"),
c("Constipation","High Normal"),
c("Constipation","Diarrhea"),
c("Low Normal","High Normal"),
c("Low Normal","Diarrhea"),
c("High Normal","Diarrhea"))
dat$bowel <- factor(dat$bowel, levels = c("Constipation", "Low Normal", "High Normal", "Diarrhea"), labels = c("Constipation", "Low Normal", "High Normal", "Diarrhea"))
#Begin preparing gut microbiome data results for plotting
#Import gut df for preprocessing:
df <- read.csv("gut_full.csv", check.names=F)
colnames(df) = gsub("nan", "Unclassified", colnames(df))
comparisons = list(c("Constipation","High Normal"),c("Low Normal","High Normal"),c("Diarrhea","High Normal"))
df$bowel <- factor(df$bowel, levels = c(1,2,3,4), labels = c("Constipation", "Low Normal", "High Normal", "Diarrhea"))
df_otus <- dplyr::select(df, -c("public_client_id","bowel","vendor_dashboard","sex","age","BMI_CALC","eGFR"))
df_otus <- as.data.frame(clr(as.matrix(df_otus)))
df_select <- dplyr::select(df, c(1:7))
df_otus <- cbind(df_select,df_otus)
df <- df_otus
dat <- df_otus
sam.n <- function(x){return(c(y = mean(x), label = length(x)))}
# ggplot2 plots
library(gplots)
Attaching package: ‘gplots’
The following object is masked from ‘package:stats’:
lowess
library(ggrepel)
library(ggbreak)
ggbreak v0.1.0
If you use ggbreak in published research, please cite the following paper:
S Xu, M Chen, T Feng, L Zhan, L Zhou, G Yu. Use ggbreak to effectively
utilize plotting space to deal with large datasets and outliers. Frontiers in
Genetics. 2021, 12:774846. doi: 10.3389/fgene.2021.774846
library(sommer)
Loading required package: Matrix
Attaching package: ‘Matrix’
The following objects are masked from ‘package:tidyr’:
expand, pack, unpack
Loading required package: crayon
Attaching package: ‘crayon’
The following object is masked from ‘package:rlang’:
chr
The following object is masked from ‘package:ggplot2’:
%+%
[]==================================================================[]
[] Solving Mixed Model Equations in R (sommer) 4.1.8 (2022-09-01) []
[] ------------- Multivariate Linear Mixed Models -------------- []
[] Author: Giovanny Covarrubias-Pazaran []
[] Published: PLoS ONE 2016, 11(6):1-15 []
[] Dedicated to the University of Chapingo and UW-Madison []
[] Type 'vignette('v1.sommer.quick.start')' for a short tutorial []
[] Type 'citation('sommer')' to know how to cite sommer []
[]==================================================================[]
sommer is updated on CRAN every 4-months due to CRAN policies
Newest source is available at https://github.com/covaruber/sommer
To install type: library(devtools); install_github('covaruber/sommer')
# Create the dfs to store BMF metadata and p values for each hit
const <- c()
const_p <- c()
low <- c()
low_p <- c()
diarrhea <- c()
diarrhea_p <- c()
taxa_names <- c()
family_names <- c()
genus_names <- c()
likelihood <- c()
adj_p <- c()
names(dv_analysis$p) = gsub("nan", "Unclassified", names(dv_analysis$p))
#Add all the p values and metadata to the right dfs:
for (i in 1:length(dv_analysis$all_models[names(dv_analysis$p)])) {
const_p[i] <- dv_analysis$all_models[[i]]$coefficients[2,4]
const[i] <- dv_analysis$all_models[[i]]$coefficients[2,1]
low_p[i] <- dv_analysis$all_models[[i]]$coefficients[3,4]
low[i] <- dv_analysis$all_models[[i]]$coefficients[3,1]
diarrhea_p[i] <- dv_analysis$all_models[[i]]$coefficients[4,4]
diarrhea[i] <- dv_analysis$all_models[[i]]$coefficients[4,1]
taxa_names[i] <- names(dv_analysis$p)[i]
family_names[i] <- strsplit(taxa_names[i],'.',fixed=TRUE)[[1]][5]
genus_names[i] <- strsplit(taxa_names[i],'.',fixed=TRUE)[[1]][6]
likelihood[i] <- dv_analysis$all_models[[i]]$logL
adj_p[i] <- dv_analysis$p_fdr[[i]]
}
genus_names[i] <- strsplit(taxa_names[i],'.',fixed=TRUE)[[1]][6]
#Create final p-value df:
p_df <- bind_cols(taxa_names,family_names,genus_names, likelihood, adj_p, const,const_p,low,low_p,diarrhea,diarrhea_p)
New names:
• `` -> `...1`
• `` -> `...2`
• `` -> `...3`
• `` -> `...4`
• `` -> `...5`
• `` -> `...6`
• `` -> `...7`
• `` -> `...8`
• `` -> `...9`
• `` -> `...10`
• `` -> `...11`
names(p_df) <- c("Genera","Family","Genus","LogL","Adj.P","Const.Beta","Const.P.Val","Low.Beta","Low.P.Val","Diarrhea.Beta","Diarrhea.P.Val")
p_df$Combined <- paste(p_df$Family,p_df$Genus)
ab <- colSums(df_otus[names(df_otus) %in% paste("taxa_",names(dv_analysis$p),sep="")])/colSums(!!df_otus[names(df_otus) %in% paste("taxa_",names(dv_analysis$p),sep="")])
p_df$Mean.Abundance <- ab[paste("taxa_",p_df$Genera,sep="")]
p_df$Const.Adj.P.Val <- p.adjust(const_p, method = "fdr", n = length(const_p))
p_df$Low.Adj.P.Val <- p.adjust(low_p, method = "fdr", n = length(low_p))
p_df$Diarrhea.Adj.P.Val <- p.adjust(diarrhea_p, method = "fdr", n = length(diarrhea_p))
p_df <- p_df[order(-p_df$Mean.Abundance),]
set <- subset(p_df[which(p_df$Adj.P < 0.05),], Genus != 'Unclassified')
set <- subset(set[order(-set$Mean.Abundance),], Genus != 'Unclassified')
abundant_list <- list()
abundant_genus <- list()
abundant <- paste("taxa_",set[order(-set$Mean.Abundance),]$Genera,sep="")
for (i in 1:68) {
abundant_list[i] <- paste(strsplit(abundant[i],'.',fixed=TRUE)[[1]][5],strsplit(abundant[i],'.',fixed=TRUE)[[1]][6])
abundant_genus[i] <- paste(strsplit(abundant[i],'.',fixed=TRUE)[[1]][6])
}
#top 10 most abundant non-unclassified hits
abundant_trunc <- set[match(abundant_list[!is.na(set$Combined[match(abundant_list,set$Combined)])],set$Combined),][which(set$Adj.P<0.05),]$Combined[1:10]
`%!in%` <- Negate(`%in%`)
t <- set[set$Combined %!in% abundant_trunc,]
#create the set of hits that are top 10 most abundant, Akkermansia, and top 9 most significant = 20 hits:
span <- rbind(set[order(set[-order(set[!set$Combined %in% abundant_trunc,]$Adj.P),][which(set$Adj.P<0.05),]$Mean.Abundance %!in% abundant_trunc),][1:10,],
t[match('Akkermansiaceae Akkermansia',t$Combined),],
t[order(t$Const.Adj.P.Val),][1:9,])
span$Letters <- LETTERS[1:20]
p_df$Letters <- c(NA)
test <- rbind(span,p_df[p_df$Combined %!in% span$Combined,])
test <- test[!duplicated(test$Genera),]
p_df <- test
p_df <- p_df[order(p_df$Adj.P),]
p_df <- p_df[order(-p_df$Mean.Abundance),]
#Pielou's Evenness:
model = y ~ x
df_num <- merge(richness, df_otus, by="public_client_id")
df_num$bowel <- factor(df_num$bowel, order = TRUE)
d1_plot <- lm(Pielou ~ bowel, data = df_num)
pval1 = summary(d1_plot)$coefficients[2,4]
r1 = summary(d1_plot)$adj.r.squared
pval1
[1] 0.01814568
r1
[1] 0.01204487
d1 <- ggplot(data = df_num, aes(x = factor(bowel, level = c("Constipation", "Low Normal", "High Normal", "Diarrhea"), labels = c("Constipation", "Low Normal", "High Normal", "Diarrhea")), y = Pielou, group = factor(bowel))) +
scale_x_discrete(bquote(atop(~italic("adj R")^2~" = "~.(formatC(as.numeric(r1),3)),~italic("P value")~" = "~.(formatC(as.numeric(pval1),3)))))+
geom_beeswarm(aes(color = bowel), cex = 0.5) +
geom_boxplot(alpha=0) +
ylim(0.005,0.05)+
geom_smooth(method = "lm", formula = y ~ x, aes(group = 1)) +
ggtitle(label = "C)\nPielou's Evenness vs BMF") +
xlab("Bowel Movement Frequency (BMF)") +
ylab("Pielou's Evenness") +
scale_colour_discrete(name="BMF Category", limits = c("Constipation", "Low Normal", "High Normal", "Diarrhea"), labels=c("Constipation", "Low Normal", "High Normal", "Diarrhea")) +
guides(colour = guide_legend(override.aes = list(size=7))) +
theme(plot.title = element_text(size=10), legend.text = element_text(size = 10), legend.title = element_text(size = 10), axis.text.x = element_blank(), axis.title.y = element_text(size = 10))
d1
Warning: Removed 1 rows containing non-finite values (stat_boxplot).
Warning: Removed 1 rows containing non-finite values (stat_smooth).
Warning: Removed 1 rows containing missing values (position_beeswarm).
#Shannon Diversity:
model = y ~ x
d2_plot <- lm(Shannon ~ bowel, data = df_num)
pval2 = summary(d2_plot)$coefficients[2,4]
r2 = summary(d2_plot)$adj.r.squared
pval2
[1] 0.005892245
r2
[1] 0.01384
d2 <- ggplot(data = df_num, aes(x = factor(bowel, level = c("Constipation", "Low Normal", "High Normal", "Diarrhea"), labels = c("Constipation", "Low Normal", "High Normal", "Diarrhea")), y = Shannon, group = bowel)) +
scale_x_discrete(bquote(atop(~italic("adj R")^2~" = "~.(formatC(as.numeric(r2),3)),~italic("P value")~" = "~.(formatC(as.numeric(pval2),3)))))+
geom_beeswarm(aes(color = factor(bowel)), cex = 0.5) +
geom_boxplot(alpha=0) +
geom_smooth(method = "lm", formula = y ~ x, aes(group = 1)) +
ggtitle(label = "B)\nShannon Diversity vs BMF") +
xlab("Bowel Movement Frequency (BMF)") +
ylab("Shannon Diversity") +
scale_colour_discrete(name="BMF Category", limits = c("Constipation", "Low Normal", "High Normal", "Diarrhea"), labels=c("Constipation", "Low Normal", "High Normal", "Diarrhea")) +
guides(colour = guide_legend(override.aes = list(size=7))) +
theme(plot.title = element_text(size=10), legend.text = element_text(size = 10), legend.title = element_text(size = 10), axis.text.x = element_blank(), axis.title.y = element_text(size = 10))
d2
#Observed ASVs:
model = y ~ x
d3_plot <- lm(Observed ~ bowel, data = df_num)
pval3 = summary(d3_plot)$coefficients[2,4]
r3 = summary(d3_plot)$adj.r.squared
pval3
[1] 0.0009015642
r3
[1] 0.0208449
d3 <- ggplot(data = df_num, aes(x = factor(bowel, level = c("Constipation", "Low Normal", "High Normal", "Diarrhea"), labels = c("Constipation", "Low Normal", "High Normal", "Diarrhea")), y = Observed, group = bowel))+
scale_x_discrete(bquote(atop(~italic("adj R")^2~" = "~.(formatC(as.numeric(r3),3)),~italic("P value")~" = "~.(formatC(as.numeric(pval3),3)))))+
geom_beeswarm(aes(color = factor(bowel)), cex = 0.5) +
geom_boxplot(alpha=0) +
geom_smooth(method = "lm", formula = y ~ x, aes(group = 1)) +
ggtitle(label = "A)\nObserved ASVs vs BMF") +
xlab("Bowel Movement Frequency (BMF)") +
ylab("Observed ASVs") +
scale_colour_discrete(name="BMF Category", limits = c("Constipation", "Low Normal", "High Normal", "Diarrhea"), labels=c("Constipation", "Low Normal", "High Normal", "Diarrhea")) +
guides(colour = guide_legend(override.aes = list(size=7))) +
theme(plot.title = element_text(size=10), legend.text = element_text(size = 10), legend.title = element_text(size = 10), axis.text.x = element_blank(), axis.title.y = element_text(size = 10))
d3
PSO <- ggarrange(d3, d2, d1, legend = "top", align = "hv", common.legend = TRUE, widths = c(7,7,7), heights = c(20,20,20), nrow = 1, ncol = 3)+ theme(plot.margin = margin(0,0,0,0))
Warning: Removed 1 rows containing non-finite values (stat_boxplot).
Warning: Removed 1 rows containing non-finite values (stat_smooth).
Warning: Removed 1 rows containing missing values (position_beeswarm).
PSO
ggsave(
"PSOvBMF.png",
plot = PSO,
device = NULL,
path = NULL,
scale = 1,
width = NA,
height = NA,
units = c("in", "cm", "mm", "px"),
dpi = 300,
limitsize = TRUE,
bg = NULL
)
Saving 6.86 x 4.24 in image
#Annotation function:
sig = function(x){
if(x < 0.001){"***"}
else if(x < 0.01){"**"}
else if(x < 0.05){"*"}
else{NA}}
#use gut df to construct column-specific df for plotting function
df <- read.csv("gut_full.csv", check.names=F)
colnames(df) = gsub("nan", "Unclassified", colnames(df))
comparisons = list(c("Constipation","High Normal"),c("Low Normal","High Normal"),c("Diarrhea","High Normal"))
df$bowel <- factor(df$bowel, levels = c(1,2,3,4), labels = c("Constipation", "Low Normal", "High Normal", "Diarrhea"))
df_otus <- dplyr::select(df, -c("public_client_id","bowel","vendor_dashboard","sex","age","BMI_CALC","eGFR"))
#refresh dataframes:
df_otus <- as.data.frame(clr(as.matrix(df_otus)))
df_select <- dplyr::select(df, c(1:7))
df_otus <- cbind(df_select,df_otus)
df <- df_otus
dat <- df_otus
#Initialize vectors:
df$p.value <- NA
df$pval <- NULL
df$p.val <- NULL
counter <<- 1
#Plotting:
# Special Graphing Function:
########################################
test = function(x,y,z,j) {
work <- NULL
results <- NULL
x = NULL
y = NULL
if (counter == 1) {z = comparisons[[1]][1]}
else if (counter == 2) {z = comparisons[[2]][1]}
else if (counter > 2) {z = comparisons[[3]][1]}
temp <- data.frame("name" = NA, "p.value" = NA, "bowel" = NA)
for (i in 1:20) {
temp$name <- paste("genus_",span$Genera[[i]],sep="")
temp$p.value <- span$Const.Adj.P.Val[[i]]
temp$bowel <- "Constipation"
work <- rbind(work,temp)
temp$p.value <- span$Low.Adj.P.Val[[i]]
temp$bowel <- "Low Normal"
work <- rbind(work,temp)
temp$p.value <- span$Diarrhea.Adj.P.Val[[i]]
temp$bowel <- "Diarrhea"
work <- rbind(work,temp)
}
if (z[[1]][1] == comparisons[[1]][1]) {results <- list(p.value = work$p.value[3*(j-1)+1])}
else if (z[[1]][1] == comparisons[[2]][1]) {results <- list(p.value = work$p.value[3*(j-1)+1+1]) }
else if (z[[1]][1] == comparisons[[3]][1]) {results <- list(p.value = work$p.value[3*(j-1)+1+2])}
counter<<-counter+1
if (counter > 3) {counter<<-1}
return(results)
}
########################################
#Store plots in list of plots for looping to plot next:
myplots <- list() # new empty list
df_test <- df
df_test[df_test==0] <- NA
for (genera in 1:20) {
y_name <- paste("taxa_",span$Genera[genera],sep="")
myplots[[genera]] <- local({
plotlim_lower = min(
df_test[!is.na(df_test[y_name]),][y_name])
plotlim_upper = max(
df_test[!is.na(df_test[y_name]),][y_name])
plotlim_bar = plotlim_lower - 4
plotlim_margin = 6
genera <- genera
plt <- ggplot(data = df_test, aes(x = bowel, y = .data[[y_name]], group = bowel)) +
scale_x_discrete(guide = guide_axis(n.dodge = 2))+
geom_jitter(aes(color = bowel), size = 0.1, cex = 0.05) +
geom_boxplot(data = df_test, alpha=0.0,outlier.shape = NA) +
theme(text = element_text(size = 9)) +
ggtitle(label = paste(span$Family[genera],"\n",span$Genus[genera],sep="")) +
geom_signif(data = df_test, comparisons = comparisons, map_signif_level = sig, test = 'test', test.args = list(z = comparisons, j = genera), y_position = plotlim_bar,step_increase = 0.15, size = 0.5,
textsize = 1.5,
tip_length = c(0,0)) +
coord_cartesian(ylim=c(plotlim_lower,plotlim_upper),clip="off")+
labs(color = "BMF Category", y = ifelse((genera == 1 | genera == 6 | genera == 11 | genera == 16),"CLR Abundance","")) +
guides(colour = guide_legend(override.aes = list(size=7), title.position = 'left', nrow = 1, ncol = 4)) +
theme(plot.margin = unit(c(0,0,plotlim_margin,0), "cm"),
plot.title = element_text(size=5.75),
legend.title = element_text(size=10),
plot.subtitle = element_text(size=10),
legend.text = element_text(size=7),
axis.text.x = element_blank(),
axis.text.y = element_text(size=7),
axis.title.y = element_text(size=7),
axis.title.x = element_blank(),
aspect.ratio = 0.95)
})
}
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
#Truncate to the top 20 including Akkermansia, most abundant, and most significant:
title <- list()
for(i in seq(1:20)) { title[i] <- myplots[[c(i)]]$labels$title }
plots1 <- myplots[title %in% paste(
span$Family,
"\n",
span$Genus,sep="")][1:20]
plots_final <- plots1[!sapply(plots1,is.null)]
#Arrange the plots:
final1 <- ggarrange(plotlist = plots_final[1:5], labels = LETTERS[1:5], legend = "top", align = "hv", font.label = list(size = 8), common.legend = TRUE, nrow = 1, ncol = 5) + theme(plot.margin = unit(c(0,0,0,0), "cm"))
Warning: Removed 3 rows containing non-finite values (stat_boxplot).
Warning: Removed 3 rows containing non-finite values (stat_signif).
Warning: Removed 3 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 3 rows containing non-finite values (stat_boxplot).
Warning: Removed 3 rows containing non-finite values (stat_signif).
Warning: Removed 3 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 56 rows containing non-finite values (stat_boxplot).
Warning: Removed 56 rows containing non-finite values (stat_signif).
Warning: Removed 56 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 2 rows containing non-finite values (stat_boxplot).
Warning: Removed 2 rows containing non-finite values (stat_signif).
Warning: Removed 2 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 85 rows containing non-finite values (stat_boxplot).
Warning: Removed 85 rows containing non-finite values (stat_signif).
Warning: Removed 85 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 171 rows containing non-finite values (stat_boxplot).
Warning: Removed 171 rows containing non-finite values (stat_signif).
Warning: Removed 171 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
final1
final2 <- ggarrange(plotlist = plots_final[6:10], labels = LETTERS[6:10], legend = "top", align = "hv", font.label = list(size = 8), common.legend = TRUE, nrow = 1, ncol = 5) + theme(plot.margin = unit(c(0,0,0,0), "cm"))
Warning: Removed 113 rows containing non-finite values (stat_boxplot).
Warning: Removed 113 rows containing non-finite values (stat_signif).
Warning: Removed 113 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 113 rows containing non-finite values (stat_boxplot).
Warning: Removed 113 rows containing non-finite values (stat_signif).
Warning: Removed 113 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 98 rows containing non-finite values (stat_boxplot).
Warning: Removed 98 rows containing non-finite values (stat_signif).
Warning: Removed 98 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 607 rows containing non-finite values (stat_boxplot).
Warning: Removed 607 rows containing non-finite values (stat_signif).
Warning: Removed 607 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 171 rows containing non-finite values (stat_boxplot).
Warning: Removed 171 rows containing non-finite values (stat_signif).
Warning: Removed 171 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 4 rows containing non-finite values (stat_boxplot).
Warning: Removed 4 rows containing non-finite values (stat_signif).
Warning: Removed 4 rows containing missing values (geom_point).
Warning: Removed 7 rows containing missing values (geom_signif).
final2
final3 <- ggarrange(plotlist = plots_final[11:15], labels = LETTERS[11:15], legend = "top", align = "hv", font.label = list(size = 8), common.legend = TRUE, nrow = 1, ncol = 5) + theme(plot.margin = unit(c(0,0,0,0), "cm"))
Warning: Removed 762 rows containing non-finite values (stat_boxplot).
Warning: Removed 762 rows containing non-finite values (stat_signif).
Warning: Removed 762 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 762 rows containing non-finite values (stat_boxplot).
Warning: Removed 762 rows containing non-finite values (stat_signif).
Warning: Removed 762 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 271 rows containing non-finite values (stat_boxplot).
Warning: Removed 271 rows containing non-finite values (stat_signif).
Warning: Removed 271 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 591 rows containing non-finite values (stat_boxplot).
Warning: Removed 591 rows containing non-finite values (stat_signif).
Warning: Removed 591 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 147 rows containing non-finite values (stat_boxplot).
Warning: Removed 147 rows containing non-finite values (stat_signif).
Warning: Removed 147 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 240 rows containing non-finite values (stat_boxplot).
Warning: Removed 240 rows containing non-finite values (stat_signif).
Warning: Removed 240 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
final3
final4 <- ggarrange(plotlist = plots_final[16:20], labels = LETTERS[16:20], legend = "top", align = "hv", font.label = list(size = 8), common.legend = TRUE, nrow = 1, ncol = 5) + theme(plot.margin = unit(c(0,0,0,0), "cm"))
Warning: Removed 592 rows containing non-finite values (stat_boxplot).
Warning: Removed 592 rows containing non-finite values (stat_signif).
Warning: Removed 592 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 592 rows containing non-finite values (stat_boxplot).
Warning: Removed 592 rows containing non-finite values (stat_signif).
Warning: Removed 592 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 40 rows containing non-finite values (stat_boxplot).
Warning: Removed 40 rows containing non-finite values (stat_signif).
Warning: Removed 40 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 331 rows containing non-finite values (stat_boxplot).
Warning: Removed 331 rows containing non-finite values (stat_signif).
Warning: Removed 331 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 680 rows containing non-finite values (stat_boxplot).
Warning: Removed 680 rows containing non-finite values (stat_signif).
Warning: Removed 680 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 20 rows containing non-finite values (stat_boxplot).
Warning: Removed 20 rows containing non-finite values (stat_signif).
Warning: Removed 20 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
final4
final_main <- ggarrange(plotlist = plots_final[1:20], labels = LETTERS[1:20], legend = "top", align = "hv", common.legend = TRUE, nrow = 4, ncol = 5)
Warning: Removed 3 rows containing non-finite values (stat_boxplot).
Warning: Removed 3 rows containing non-finite values (stat_signif).
Warning: Removed 3 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 3 rows containing non-finite values (stat_boxplot).
Warning: Removed 3 rows containing non-finite values (stat_signif).
Warning: Removed 3 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 56 rows containing non-finite values (stat_boxplot).
Warning: Removed 56 rows containing non-finite values (stat_signif).
Warning: Removed 56 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 2 rows containing non-finite values (stat_boxplot).
Warning: Removed 2 rows containing non-finite values (stat_signif).
Warning: Removed 2 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 85 rows containing non-finite values (stat_boxplot).
Warning: Removed 85 rows containing non-finite values (stat_signif).
Warning: Removed 85 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 171 rows containing non-finite values (stat_boxplot).
Warning: Removed 171 rows containing non-finite values (stat_signif).
Warning: Removed 171 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 113 rows containing non-finite values (stat_boxplot).
Warning: Removed 113 rows containing non-finite values (stat_signif).
Warning: Removed 113 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 98 rows containing non-finite values (stat_boxplot).
Warning: Removed 98 rows containing non-finite values (stat_signif).
Warning: Removed 98 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 607 rows containing non-finite values (stat_boxplot).
Warning: Removed 607 rows containing non-finite values (stat_signif).
Warning: Removed 607 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 171 rows containing non-finite values (stat_boxplot).
Warning: Removed 171 rows containing non-finite values (stat_signif).
Warning: Removed 171 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 4 rows containing non-finite values (stat_boxplot).
Warning: Removed 4 rows containing non-finite values (stat_signif).
Warning: Removed 4 rows containing missing values (geom_point).
Warning: Removed 7 rows containing missing values (geom_signif).
Warning: Removed 762 rows containing non-finite values (stat_boxplot).
Warning: Removed 762 rows containing non-finite values (stat_signif).
Warning: Removed 762 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 271 rows containing non-finite values (stat_boxplot).
Warning: Removed 271 rows containing non-finite values (stat_signif).
Warning: Removed 271 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 591 rows containing non-finite values (stat_boxplot).
Warning: Removed 591 rows containing non-finite values (stat_signif).
Warning: Removed 591 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 147 rows containing non-finite values (stat_boxplot).
Warning: Removed 147 rows containing non-finite values (stat_signif).
Warning: Removed 147 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 240 rows containing non-finite values (stat_boxplot).
Warning: Removed 240 rows containing non-finite values (stat_signif).
Warning: Removed 240 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 592 rows containing non-finite values (stat_boxplot).
Warning: Removed 592 rows containing non-finite values (stat_signif).
Warning: Removed 592 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 40 rows containing non-finite values (stat_boxplot).
Warning: Removed 40 rows containing non-finite values (stat_signif).
Warning: Removed 40 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 331 rows containing non-finite values (stat_boxplot).
Warning: Removed 331 rows containing non-finite values (stat_signif).
Warning: Removed 331 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 680 rows containing non-finite values (stat_boxplot).
Warning: Removed 680 rows containing non-finite values (stat_signif).
Warning: Removed 680 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 20 rows containing non-finite values (stat_boxplot).
Warning: Removed 20 rows containing non-finite values (stat_signif).
Warning: Removed 20 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
ggsave(
"BMFvsGenera1.png",
plot = final1,
device = NULL,
path = NULL,
scale = 1.5,
width = NA,
height = NA,
units = c("in", "cm", "mm", "px"),
dpi = 300,
limitsize = TRUE,
bg = NULL
)
Saving 10.3 x 6.36 in image
ggsave(
"BMFvsGenera2.png",
plot = final2,
device = NULL,
path = NULL,
scale = 1.5,
width = NA,
height = NA,
units = c("in", "cm", "mm", "px"),
dpi = 300,
limitsize = TRUE,
bg = NULL
)
Saving 10.3 x 6.36 in image
ggsave(
"BMFvsGenera3.png",
plot = final3,
device = NULL,
path = NULL,
scale = 1.5,
width = NA,
height = NA,
units = c("in", "cm", "mm", "px"),
dpi = 300,
limitsize = TRUE,
bg = NULL
)
Saving 10.3 x 6.36 in image
ggsave(
"BMFvsGenera4.png",
plot = final4,
device = NULL,
path = NULL,
scale = 1.5,
width = NA,
height = NA,
units = c("in", "cm", "mm", "px"),
dpi = 300,
limitsize = TRUE,
bg = NULL
)
Saving 10.3 x 6.36 in image